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Abstract 

The effects of demographic stochasticity in the long term behaviour of endemic infectious diseases 
have been considered for long as a necessary addition to an underlying deterministic theory. The 
latter would explain the regular behaviour of recurrent epidemics, and the former the superimposed 
noise of observed incidence patterns. Recently, a stochastic theory based on a mechanism of 
resonance with internal noise has shifted the role of stochasticity closer to the center stage, by 
showing that the major dynamic patterns found in the incidence data can be explained as resonant 
fluctuations, whose behaviour is largely independent of the amplitude of seasonal forcing, and 
by contrast very sensitive to the basic epidemiological parameters. Here we elaborate on that 
approach, by adding an ingredient which is missing in standard epidemic models, the 'mixing 
network' through which infection may propagate. We find that spatial correlations have a major 
effect in the enhancement of the amplitude and the coherence of the resonant stochastic fluctuations, 
providing the ordered patterns of recurrent epidemics, whose period may differ significantly from 
that of the small oscillations around the deterministic equilibrium. We also show that the inclusion 
of a more realistic, time correlated, recovery profile instead of exponentially distributed infectious 
periods may, even in the random- mixing limit, contribute to the same effect. 
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I. INTRODUCTION 



The incidence patterns of childhood diseases in the twentieth century have been a chal- 
lenge and a preferred testing ground for epidemiological models. During the last decade, 
more sophisticated approaches building on the traditional SIR and SEIR models |l| (see 
Section II for a description) have brought considerable advances in understanding and se- 
lectirig some of the fundamental ingredients of the complex dynamics of infectious diseases 



, [^). The interplay between the system's nonlinearity and the periodic perturbation 
in seasonally forced SIR and SEIR models was shown to be a source of complex dynamics 

n n Q 

compatible with the diversity of observed incidence records ([5[, [6|, [7|). 

This body of work belongs to an essentially deterministic framework, where demographic 
stochasticity plays a secondary role, except when addressing stochastic extinction 
In this framework, the role of stochasticity is that of sustaining small amplitude oscillations 
around the deterministic system's equilibrium that follow the natural frequency given by the 
local linear approximation ([gI, [l^, [ill), or else that of promoting the switching between 
different competing attractors of the underlying deterministic model ( jl2|. 



Recently, a stochastic theory developed for a predator-prey competition model and 
then appHed to the weU-^xed SIR .node. Q has shown that the fluetuat.on power spectrum 
of the incidence time series is essentially determined, both in the presence and in the absence 
of seasonal forcing, by the resonance with internal noise of the system's natural frequency 
in the deterministic (infinite population size) limit. The cross correlation structure of this 
internal noise can be computed analytically from the transition probabilities, and it may 
shift the resonance peak away from the system's natural frequency. The amplitude and 
the coherence of these resonant fluctuations are both large for the parameter values that 
correspond to some childhood diseases, so that they may be comparable even in large systems 
to the oscillations induced by seasonal forcing. 

Apart from stochasticity, another missing ingredient of the standard epidemic models 
that has been attracting increasing attention is the host population contact structure, or 

n n 

'mixing network' [16[, [17[. Most of the research connecting networks and epidemiology 
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dates from the last five years, when many fundamental results of network theory became 

n n □ 

widely known while new ones have been derived [18| , [19[ , [20| . These ideas originated in the 
mathematics and physics communities were applied in an epidemiological setting from the 

beginning J21J , 22|], 23| prompting the interest of epidemiologists and several theoretical 

24|, 25], Q] and field epidemiology 28], 29] results. 

Here we address the effect of relaxing the random mixing assumption on the behaviour 
of the resonant fluctuations of a stochastic SIR model. We assume that the host popula- 
tion contact network may be represented by a 'small-world' network of the type introduced 
by Watts and Strogatz nearly a decade ago 211]. Analysis of real networks of social con- 
tacts SOj indicates that this is a reasonable first assumption to model the contact structure 
relevant for the propagation of airborne infections. We find that spatial correlations have 
a major effect in the enhancement of the amplitude and of the coherence of the resonant 
stochastic fluctuations, providing ordered patterns of recurrent epidemics. The shift of the 
resonance peak away from the system's natural frequency also increases significantly due 
to correlations. The enhanced amplitude and coherence of the fluctuations implies that in 
populations where a disease spreads predominantly through local infectious contacts, large 
epidemic outbursts with well defined recurrence times are generated by internal noise alone, 
without the presence of seasonal forcing. This effect is illustrated in Section II. B with a 
particular example, and studied systematically in Section III.B. 

Another assumption of the SIR model which has been challenged 311], 32|], 33| is that of 
considering a constant recovery probability during the infection period. We study the effect 
on the resonant fluctuations of switching from the (uncorrelated) exponentially distributed 
infection periods used in 15| to the opposite limit of (time correlated) constant infection 



periods. Again, we flnd that the spectrum of stochastic fluctuations found in [15| changes 
under more realistic assumptions on the recovery proflle in such a way that the oscillations 
become larger and more sharply deflned. In this case, the peak frequency shift with respect 
to the deterministic prediction is larger, and towards larger fequencies. 



Qualitatively, the effects of spatial and time correlations on the fluctuation spectrum are 
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similar: both the spatial correlations, introduced through the host population contact net- 
work, and the time correlations, introduced via constant infection periods, lead to enhanced 
more sharply defined dominant peaks. There arc analogues of this behaviour in many sys- 
tems studied in statistical physics. The peak frequency shift, however, is more pronounced 
and has a different sign for the model with constant infection period. As discussed in Sec- 
tions II. B and III.C, this quahtative difference can be understood on the basis of an effective 
deterministic description. 

These results show that the classical understanding of the incidence time patterns of 
endemic infectious diseases, which is mainly based on a seasonally forced deterministic de- 
scription, is clearly insufficient to have a correct description of such fiuctuations. These 
may be purely stochastic, and the diversity of incidence patterns found in real data for the 
same disease in different populations can be understood in this framework as an effect of 
population size and contact structure. 

Moreover, the approximate period of the recurrent epidemics driven by internal noise in 
the presence of spatial and/or time correlations differs significantly from the period computed 
in the usual deterministic approach. This is an important cautionary note with regard 
to the use, in the framework of the deterministic description, of the recurrence period to 
help assess estimates of epidemiological parameters. The breakdown of the assumptions 
of random mixing of the population and/or of constant recovery rate during the infectious 
period implies important corrections to the dominant frequency of the fluctuations. 

II. METHODS 

A. Susceptible-Infectious-Recovered (SIR) dynamics in a randomly mixed dis- 
crete population 

One of the simplest epidemiological models one can consider is based on dividing the 
whole population in three classes of individuals, the susceptible, the infectious, and the 
permanently recovered. It is suitable to study the infection dynamics of diseases that confer 
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long lasting immunity, such as childhood infections, provided that we take into account in 
the long term dynamics the replenishment of susceptibles in the population through births. 

The deteministic description of the basic Susceptible-Infectious-Recovered (SIR) model 
in a closed population where renewal of susceptibles occurs through births with a constant 
birth rate and death rate /i is given by the equations 

s = yu(l — s) — Psi 

i = /3si - (7 + /i)^ (1) 

where s, i are the densities of susceptible and infectious individuals, respectively, 7 the 
recovery rate of the disease and (3 its transmissibility. A crucial parameter for the behaviour 
of this model is the so called basic reproductive rate, Rq = P/{'y + fi). For _Ro < 1 the 
disease dies out, while for Rq > 1, the disease is endemic in the population. In this case, 
system (JT]) has a non trivial endemic equilibrium at s* = (7 + yu)//3, ^* = (/i//5)(l/s* — 1), 
and the small oscillations around this equilibrium have damping factor — /i/ (2s*) and period 
27r/(y^/i/5(l — s*) — (/i/2s*)2). The Susceptible-Exposed-Infectious-Recovered (SEIR) model 
considers an additional class of infected, but not yet infectious, individuals. 

In a randomly mixed discrete population of individuals, SIR dynamics may be de- 
scribed as a continuous time Markov process on a population divided into three classes, 
susceptibles, infectious and recovered. The state of the system is characterized by the num- 
bers S, I and R = N — S~I, of individuals in each of the three classes, and the events of 
infection, death, birth and recovery correspond to the following transitions and transition 
rates 

T[(5, /) ^ (S - 1, / + 1)] = (3SI/N infection 

T[{S, /) ^ (5 — 1, /)] = fiS death of a susceptible 
T[{S, I) ^ {S, I — 1)] = fil death of an infectious 

T[{S, I) {S, I)] = fi{N — S — I) death of a recovered 
T[{S, I)^{S + 1, /)] = 12N birth 

T[(5,/) (5,/- 1)] = 7/ recovery, (2) 
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where {S, I) (S", /') denotes the transition from state {S, I) to state (S", /') and T[{S, I) —>■ 
{S', I')] the corresponding transition rate. 

As detailed in ISj, a good approximation for the power spectrum of the stochastic fluc- 
tuations around the stationary population numbers can be computed analytically from the 
linear Fokker-Planck equation obtained from the next to leading order terms in van Kam- 
pen's expansion of the master equation associated to The leading order terms of the 
expansion yield the deterministic equations ([T]), that describe the behaviour of the system 
in the limit of infinite populations. 

Following this approach, we have computed the power densities Ps and Pi of the fluctu- 
ations of susceptibles and infectious individuals for process ([2]), scaled by the square root of 
the system size V^, as a function of the angular frequency ui, 



Psiu) = 2/i (^1 
Piioo) = 2/i (^1 

where D = — j — fi) and T = — /9yu/(7 + fi) = —Rofi, denoting as before by Rq the basic 
reproductive rate of the disease. The parameters D and T are equal to the determinant and 
to the trace, respectively, of the linear approximation of ([H) at the endemic equilibrium, and 
they have a simple dynamical interpretation: D is the square of the frequency of the small 
oscillations around the equilibrium when the damping is small, and T is twice the damping 
factor of these oscillations. 

Equations ([S]) are independent of because the dependence of the amplitude of the 
fluctuations on system size has been scaled out, and they describe exactly the stochastic 
process ([2]) in the limit of large A^. These power spectra are resonant like, showing that the 
amplitude of the fluctuations as a function of their time scale is governed by a mechanism of 
resonance of internal noise with the system's natural frequency -\/zJ- The resonance peak will 
be shifted away from \/D, depending on the values of T and on the terms in the numerator of 
(E]), which are determined by the cross correlation structure of the internal noise. However, 



1 + + - /3/(7 + /i) + (/?/(7 + /i))^) 
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when T is small, the shift in the resonance peak with respect to v -D is also small. 

The relevance and universality of this type of 'endogenous' stochastic resonance for eco- 
log,cal systems i„ gene,^ was first argued irr Q. The complete analyt.c description of the 
phenomenon given in [ij] , [15| relies on the random mixing assumption and on the constant 
recovery rate assumption. The effect of relaxing these assumptions, crucial in more realistic 
settings, must be tested through stochastic simulations. 



B. The SIR model on dynamic small- world networks 



Complex network theory [l^, [l9|], 20|] focuses on abstract network models such as "small- 
world" networks, which interpolate between regular lattices and random graphs, and scale- 
free networks, exhibiting a power law distribution of the connectivity. In most, if not all, of 
the theoretical work in contact network epidemiology, the host population contact network 
is represented by one of these models. The choice of the right idealized network type, which 
will be disease specific, depends on the availability of data on disease causing contacts which 
is hard to get and difficult to interpret 16|. At this point, an important concern of contact 
network epidemiology is to collect data and build appropriate network models for different 
transmission mechanisms 30 1. 

As a first step to understanding the role of network structure correlations on the spectrum 
of stochastic fiuctuations, we have modelled the mixing network of the populations as a 
small world network 2l|] built over a square lattice with 12 nearest neighbours per node. In 



these models a fraction of the links of the lattice is randomized by connecting nodes, with 
probability p, with any other node. These non-local connections are chosen randomly for 
each event, instead of being fixed in a frozen, partially random, link configuration. This 
version of the small world network model, which has been dubbed 'dynamic' or 'annealed' 
in the literature ( 3J], 35|], 36|]), is motivated by the nature of the occasional social contacts 



the model tries to represent. For p = 0, each node interacts only with its nearest neighbours 
on the lattice, as in ordinary representation of spatial structure. For p = 1, the network of 
interactions is a random graph, where every pair of nodes, independently of the distance on 
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the lattice between the two nodes, has the same probabihty of being connected. Random 
graphs have the property that the average path length, or average number of connections 
of the shortest path between two nodes, is 'small', i. e., of the order of the logarithm of 
the total number of nodes. For a range of p between and 1 the network exhibits 'small 
world' behaviour, where predominantly local interactions (as in lattices) coexist with a 
short average path length (as in random graphs). Analysis of real networks [30j] reveals 
the existence of small world behaviour in many interaction networks, including networks of 
social contacts. 

An important consequence of the spatial correlations introduced by predominantly local 
contacts is what is called infective screening, or infective clustering. If all the infected 
neighbours of an infected node have many neighbours in common, each of them will be 
connected to a number of susceptibles which is smaller that the average number of susceptible 
neighbours per node, and infection will be less likely than in a randomly mixed population 
with the same density of infectious and of susceptibles. 

The qualitative analysis of time series of SIR and SEIR stochastic simulations on dynamic 
small world networks for different values of the small world parameter p shows that indeed the 
infection process becomes less and less efficient as p decreases and infective screening becomes 
more pronounced 37|], 38j. The global effect of infective clustering may be quantified in 
terms of the value of Peff, defined as the average number of new infections per time step and 
infective individual divided by the density of susceptibles. For p = I, Peff coincides with 
the transmissibility (3 but, as p decreases, f3eff also decreases. Another qualitative effect of 
spatial correlations reported in [s^], jsS] is that they have a major effect on the enhancement 
of the amplitude of stochastic fluctuations, which become more and more pronounced as p 
decreases. 

In order to quantify the effect of network structure on resonant stochastic fluctuations, 
we have computed the power spectrum of SIR time series on dynamic small world networks 
for several values of p, and we have compared it with the analytic power spectrum given by 
([3]) with Peff, the effective transmissibility, instead of P (see the Appendix for the details 



8 




0.01 0.02 0.03 0.04 0.05 0.06 0.01 0.02 0.03 0.04 0.05 0.06 0.01 0.02 0.03 0.04 0.05 0.06 



Angular Frequency Angular Frequency Angular Frequency 

FIG. 1: Fluctuation power spectra of infectives time series for the SIR model on dynamic networks 
of = 10^ nodes for small world parameter p equal to (a) 1, (b) 0.4, and (c) 0.2; the networks 
are built over regular 1000 x 1000 lattices with 12 first neighbours and periodic boundary condi- 
tions. The full lines are the averaged numerical power spectra of 400 stationary time series, each 
32768 timesteps long. The dashed lines are the analytic power spectra given by (3), where for p 
smaller than 1 the value of the parameter (3 has been taken as f3cs, defined as the actual number 
of new infections per time step and infective individual divided by the density of susceptibles. 
Parameter values: /x = 6 x 10"'^, (3 = 1.32, 7 = 1/8. 

of the simulations). The results are shown in Figure 1. The full lines are the averaged 
numerical power spectra of the stationary time series, and the dashed lines are the analytic 
power spectra given by ([3]) with the correction due to the effective transmissibility. We see 
that, as p decreases, the resonant fluctuations are larger and more coeherent. Moreover, 
in the small world regime, where local correlations become important, this effect is much 
more pronounced than in the predictions of an effective (corrected for infective screening) 
randomly mixed model, represented by the dashed lines in Figure 1. However, this effective 
model does describe satisfactorily the shift to the left in the peak frequency, which means 
that this can be understood as a consequence of the reduced effective transmissibility. 

The enhanced amplitude and coherence of the fluctuations implies that a typical time 
series will exhibit noisy but regular incidence oscillations, as shown in Section III.B. These 
sustained oscillations are of a purely stochastic nature, and they disappear in the infinite 
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population limit. 

The influence of the recovery proflle on the behaviour of the system has been discussed in 
3l| . [331], [391]. In the standard SIR stochastic model, to which the analytic description ([3]) 
applies, the event of recovery occurs at a flxed rate during the infection, and the recovery time 
is exponentially distributed around the average infectious period r = I/7, or in other words 
uncorrelated. We have computed the power spectrum of the time series obtained from SIR 
simulations in randomly mixed populations where instead of uncorrelated infection periods 
the recovery proflle was taken in the opposite limit of constant infection period r or strongly 
time correlated infections. The results are shown in Figure 2. The full line is the averaged 
numerical power spectra of the stationary time series, and the dashed line is the analytic 
power spectra given by ([3]) with 7 = 1/r. We see that switching to more realistic (time 
correlated) recovery proflles leads to a similar effect of enhancing the amplitude and the 
coherence of the resonant stochastic fluctuations predicted by the theory developed in |l5 ]. 
In addition, there is a large shift to the right of the dominant frequency of the fluctuations, 
relative to the peak frequency predicted for the model with stochastic recovery [15] . Again, 
the shift in the peak frequency and its sign can be understood in terms of an effective 
deterministic description, as discussed in Section III.C. 

Qualitatively, it is to be expected that correlations in the dynamics of the system's 
components should enhance the global density fluctuations. To assess quantitatively the 
effect of spatial and temporal correlations on the resonant fluctuations spectrum in the 
modelling of infectious childhood diseases, we have performed systematic simulations in a 
region of parameter space that includes the values for measles, chicken pox, rubella, pertussis 
and mumps according to published and estimated data for the pre- vaccination period Q]. 
The amplitude and coherence of the stochastic fluctuations are measured by the overall 
amplification A, and by the coherence factor c, introduced in 15|, where analytic results 
for a randomly mixed SIR stochastic model with external infection and constant recovery 
rate were presented for the same region of parameter space. The overall amplification A is 
the integral over all the frequencies of the power density of the scaled fiuctuations, and it is 
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FIG. 2: Fluctuation power spectra of infectives time series of stochastic SIR dynamics with de- 
terministic and stochastic recovery on a randomly mixed population with N = 10^ individuals. 
Deterministic recovery occurs r timesteps after infection. The full line is the averaged numerical 
power spectra of 400 stationary time series, each 32768 timesteps long. The dashed line is the ana- 
lytic power spectra given by (3), with 7 = l/r. Parameter values: ^ = 6 x 10"*^, f3 = 1.32, r = 8. 

equal to the mean square deviation of the time series from the equilibrium values, divided 
by the system size A^. The coherence factor c is defined as the integral of the power density 
of the of the scaled fluctuations on the frequency interval that corresponds to periods within 
10% of the peak period 2TT/ujpeak, divided by A. It measures the relative contribution to 
the overall amplification of stochastic fluctuations that are distributed around the dominant 
period. 

For the randomly mixed case p = 1 with stochastic recovery, both A and c can be 
computed analytically from equations ([3]), which are exact in the limit of large A^. In the 
remaining cases there is no analytic description, and the overall amplification A is computed 
as the ensemble average, over several runs, of the integral over all sampled frequencies of 
the power density of the stationary time series of the scaled fluctuations. The coherence 
factor c is computed in a similar way on the prescribed frequency interval. The details of 
the analytical and numerical computations are given in the Appendix. 

For the model with stochastic recovery, we have also computed, for several values of the 
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small world parameter p, the peak shift factor s, defined as the distance between the actual 
peak frequency and the natural frequency of the system in the deterministic description ([T]), 
\/D, divided by \/D. It measures the relative frequency shift of the dominant frequency due 
to the various ingredients that are missing in the deterministic description (resonance with 
correlated internal noise, and contact network structure). 



III. RESULTS 



A. Randomly mixed model: finite size effects 

The randomly mixed model depends on the demographic and epidemiological parameters 

n 

fi, (3 and 7. Following [15'], we have taken fixed and equal to 5.5 x 10 , and we have 
taken the reduced variables in the parameter plane (/i/7,/?/7). These reduced parameters 
have an immediate epidemiological interpretation, /i/7 measures the acuteness, or relative 
time scale, of the disease, and /3/7 ~ -Rq- 
The values of the overall amplification 

1 /•+00 

A=- Pi{uj)duj (4) 



7r Jo 

and of the coeherence 

1 r^peak/O-9 

c = ^ / Pi{u)duj (5) 

of the infectives power spectrum Pj in are shown in Figure 3. (a) and 3.(b) for 441 points 
in parameter space. The values of A are normalized by the largest overall amplification 
0.1902. We see that, for the basic SIR model, A is essentially determined by Rq, and that 
both the overall amplification and the coherence increase as Ro decreases with 7 fixed. 
For each Rq the stochastic oscillations become more and more coherent as 7 increases. The 
symbols mark the parameter values for measles, chicken pox, rubella and pertussis according 
to different data sources for the pre- vaccination period 6|]. The numerical results for the 
overall amplification and coherence obtained from simulations of the stochastic process ([2]) 
on a population of size = 10^ are plotted in Figures 3.(c) and 3.(d) (see the Appendix for 
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FIG. 3: Overall amplification and coherence of the randomly mixed model with stochastic recov- 
ery for 21 X 21 points in parameter space. The symbols mark the parameter values for measles 
(triangles), chicken pox (circles), rubella (squares) and pertussis (diamonds) according to different 
data sources for the pre- vaccination period. Analytic results for the amplification A in (a), and 
for the coherence c in (b). The values of (a) are normalized by the largest overall amplification 
0.1902. In (c) and (d) are plotted the numerical results for the overall amplification and coherence 
obtained from simulations of the stochastic process (2) on a population of size = 10^. Shown 
in (c) are the values of the analytic plot (a) divided by the corresponding numerical values. Grey 
areas denote regions of parameter space for which less than 200 surviving runs 2^^ timesteps long 
were obtained out of 500 trials. 
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the details of the simulations and numerical computations). Each data point in Figure 3.(c) 
is the ratio of the corresponding analytic and numerical values of the overall amplification, 
while each data point in Figure 3.(d) corresponds to the numerical value of the coherence 
as defined in (I5l). The data points shown in grey correspond to parameter values where 
more than 60% of the 50 000 days long simulations ended up with zero infectives. It is clear 
that for realistic population sizes, in a large region of parameter space, there are corrections 
to the analytic results of Figure 3. (a) and 3.(b), and that the fluctuations for small (3 and 
large 7 are larger than those predicted by the analytic approach. The numerical power 
spectra obtained from stochastic simulations for = 10^ and = 50 x 10^ for values of 
(/i/7,/5/7) for which these corrections are more significant are shown in Figure 4. (a). The 
analytic power spectrum given by ([3]) overlaps the numerical curve for = 50 x 10^ within 
the resolution of the figure. A breakdown of the analytic description for A^ = 10^ can be 
detected in the overall amplification and in the peak frequency shift, as well as in the loss of 
coeherence due to the appearance of a small harmonic peak. This effect is more pronounced 
in the presence of spatial correlations, as shown in Figure 4. (b) and discussed in the following 
section. 

The limitations of the analytic description ([3]) are due to the fact that a linear theory 
based on van Kampen's method is used to approximate the master equation up to next to 
leading order tems, leaving out the contribution of terms of order 1 / y/N and higher, whose 
influence on the dynamics becomes more important as the system size and/or the stability 
of the stationary state decreases. 

B. Small world network: effects of spatial correlations 

To assess the effect of spatial correlations on the resonant fluctuations spectrum, we 
have performed, for the same 441 points in parameter space, systematic simulations of SIR 
time series on dynamic small world networks for several values of p (see the Appendix for 
the details of the simulations). The results for the overall amplification A when p = 0.6 
are shown in Figure 5. a). The values of A are normalized by the largest value of the 
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FIG. 4: (a) Numerical power spectra obtained from stochastic simulations for N = 10^ and = 
50 X 10^ (full and dotted lines, respectively) for values of {fJ-/j,P/"f) = (0.00094,7.9). The spectra 
were taken from 2^^, instead of 2^^, timesteps long simulations for increased resolution in frequency. 
The curve for N = 10^ shows a larger overall amplification and a peak frequency shift, as well as 
loss of coherence due to the appearance of a small harmonic peak, (b) Power spectra of stochastic 
simulations on a dynamic small world network with p = 0.2 and N = 10^ nodes (full line) and 
= 50 X 10^ (dotted line), for the parameter values {fi/'y,f3/^) = (0.00182,4.0) chosen in the 
region where the amplitude of the harmonic peaks is larger. Both plots are scaled to the largest 
spectral power (3960 for N = 10^ and 71000 for iV = 50 x 10^). The ratio of the heights of the 2'^'^ 
peak to the 1^* peak is approximately 7 ~ y/50 times smaller for the larger system, in line with the 
expected l/\/iV dependence. 

overall amplification, which is 0.2908. With respect to the analytic spectrum ([3]), the overall 
amplitude of the fluctuations increases by a factor of about 1.5, causing more ocasional 
extinctions. As before, the grey region corresponds to parameter values where more than 
60% of the 50 000 days long simulations ended up with zero infectives. The dependence 
of the amplitude of the fluctuations on the epidemiological parameters is qualitatively the 
same as for the randomly mixed {p = 1) numerical results of Figure 3.(c). It increases 
as -Ro decreases, and, for fixed Rq, it increases as 7 increases. Close to the extinction 
boundary, especially for low Rq we find a departure from this general trend which can be 
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FIG. 5: (a) Overall amplification of the SIR model with stochastic recovery on a dynamic small 
world network with N = 10^ individuals and p = 0.6 (see the Appendix for details of the simula- 
tions). The values of the overall amplification A are normalized in the plot by the largest amplitude 
(0.2908). Again, grey areas denote regions where only a small percentage of timeseries survive af- 
ter 2^^ timesteps (each sampled point requires the survival of at least 200 out of 500 runs), (b) 
Amplification A (diamonds) and coherence c (circles) for the point (/i/7,/?/7) = (0.00116,17.0) as 
a function of the small world parameter p. The chosen point lies approximately at the center of 
the region plotted in (a). 

due a sampling bias, as we have had to select long time series without disease extinction 
from a large ensemble of trial runs. 

The most relevant effect of spatial correlations is the increase in the amplitude of the 
fluctuations as the small world parameter p decreases. The plot of the overall amplification 
A as a function of p is shown in Figure 5.b) (data points plotted with diamonds) for parameter 
values close to those of pertussis (for other disease parameter values the behaviour of the 
amplitude of the fluctuations as a function of p is similar). There is a four-fold increase in 
the amplitude of the fluctuations for p = 0.2 with respect to the randomly mixed case p = 1. 

For the same parameter values, the plot of the coherence c as a function of p is also shown 
in Figure 5.b) (data points plotted with circles). There is an increase in the coherence of the 



16 



fluctuations as p decreases, which means that as the fluctuations get larger they also exhibit 
a more regular temporal pattern. For a fixed value of p, the coherence is uniformly very 
high in parameter space (results not shown), with changes of less than 10% in a region that 
includes all the parameter values considered for childhood infectious diseases. However, in 
a small region close to the (3/^ = 4 horizontal line and for p = 0.2 there is loss of coherence, 
due to the appearance of harmonic peaks, see Figure 4.(b), which are more pronounced than 
the ones found for p = 1, and persist for larger values of A^. Although the relative amplitude 
of the harmonic peaks seems to scale with l/\/N, the fact that they are enhanced is due to 
the presence of spatial correlations, and could be related to the existence of an oscillatory 
phase for small values of yu/7 in deterministic SIR models that include, at the simplest level, 
the effect of these correlations . Indeed, the breakdown of Van Kampen's 1 / scaling 
is expected in the oscillatory phase of this model, where the ratio of the amplitudes of the 
harmonic peaks is constant, in the infinite size limit. 

The combined effect of the spatial correlations on the amplitude and on the coherence 
of the fluctuations can be seen in Figure 6, where a typical time series for the parameter 
values of Figure 5.(b) and p = 0.2 is shown. The purely stochastic incidence peaks can be 
as high as 3500 individuals per day in a population of = 10^, with troughs of one or two 
hundred infectious, and a period of recurrence of epidemics can be clearly identified. Indeed, 
the pattern of sustained oscillations shown in Figure 6 is similar to that found in incidence 
data for measles in the prevaccination records of english cities of comparable population size 
411 1 , and it is more sharply defined than that of real time series for most childhood diseases 

n n — 

incidence data [42j, [43l |. 



44, |45|. 



Another important effect of contact network structure is the shift in the dominant fre- 
quency of the power spectrum. In Figure 7 we plot, as a function of p, the peak shift factor 
s, defined as the difference between the actual peak frequency and the natural frequency of 
the system in the deterministic description ([T]), ^/D, divided by ^/D. Again, we have chosen 
the same parameter values of Figure 5.b), (/i/7,/3/7) = (0.00116, 17.0), but the results for 
other points in parameter space are similar. We found that, as p decreases, the dominant 
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FIG. 6: Typical incidence time series for a population oi N = 10^ and small world parameter 
p = 0.2. We have taken {ii/-/, /3/-f) = (0.00116, 17.0) as in Figure 5.b). 
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FIG. 7: The peak shift factor s as a function of the small world parameter p. We have chosen 
(/x/7,/?/7) = (0.00116, 17.0), the same parameter values as in Figure 5.(b). 

frequency of the time series is shifted increasingly to the left. As discussed in Section II, this 
can be undestood as a result of the reduced effective transmissibility due to clustering of 
infectives. Since both the overall amplification and the coherence increase as p decreases, the 
time series of long term simulations will, as shown in Figure 6, exhibit recurrent epidemics 
with approximate period close to that given by the dominant frequency of the resonant 
fluctuations. However, this value can be very different from the period that corresponds to 
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FIG. 8: Overall amplification (a) and coherence (b) for the randomly mixed SIR model with 
deterministic recovery over a population of = 10^ individuals. The values of (a) are normalized 
by the largest overall amplification 0.6901. Again, grey areas denote regions of parameter space 
where less than 200 out of 500 runs survive 2^^ timesteps. 

the natural frequency of the system in the deterministic limit ([T]). 
C. Recovery profile: effects of time correlations 

The results for the overall amplification and coherence in the the randomly mixed SIR 
model with deterministic recovery are shown in Figure 8. The values of A in Figure 8. (a) 
are normalized by the largest value of the overall amplification 0.6901, which means that the 
overall amplitude of the fluctuations increases by a factor of about 3.5 with respect to the 
model with stochastic recovery. In agreement with the results of Figure 2, the enhancement 
of the peak amplitude will in general be even larger, since, as shown in Figure 8.(b), the co- 
herence of the fluctuations for the case of deterministic recovery is larger for most parameter 
values. Realistic recovery profiles will therefore be associated with time series whose fluctu- 
ations power spectra will be larger and more sharply peaked than for constant recovery rate 
models. The peak frequency shift relative to this model is towards larger frequencies, and it 
may partially cancel the effect of internal noise correlations that shifts the peak frequency 
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to the left of v-D, bringing the resulting dominant frequency closer to the natural frequency 
of the system in the deterministic description. 



In 39|, the basic SIR model ([T]) was modified to include a family, parametrized by an 
integer n, of infectious periods distributions that interpolate between the exponential distri- 
bution (for n = 0) and delta distribution associated to deterministic recovery (in the limit 
n oo). Although there is no simple expression in closed form in this case for the period 
of the damped oscillations close to the system's equilibrium, it was shown that the period 
decreases as n increases. Our findings are in agreement with the expectation of an increase 
of the dominant frequency in the limit of fixed recovery time. 

IV. DISCUSSION AND CONCLUSIONS 

The relevance of the phenomenon of resonance with internal noise for the understanding 
of the long term dynamics of childhood infections, with or without the presence of seasonal 
forcing, has been uncovered in 1^. The power spectrum of the fluctuations of the stochastic 
process of infection in a population has been given a complete analytic description, under 
the following basic assumptions: 

1. There is external infection at a small constant rate 

2. The population is randomly mixed and internal infection follows a law of mass-action 
with constant transmissibilty 

3. Recovery from the disease occurs at a constant rate 

4. The behavior of the fluctuations is well described by the lowest order van Kampen 
expansion of the master equation, 

that correspond to a stochastic open SIR model and large population sizes. 

We have extended the analytic results of [l5| to a closed stochastic SIR model, with no 
external infection, that models an isolated population, and compared these with the results 
of extensive numerical simulations. Due to stochastic extinctions, simulations of the long 
term behaviour of this model require large population sizes, and we took = 10^ as a 
typical value for the population size. We have found a good agreement between the analytic 
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description and the results of the simulations for most values of the parameters. However, 
in a region of parameter space that includes in particular measles, we have found that there 
are significant finite size corrections to the analytic power spectrum of the fiuctuations, 
which means that assumption 4 is fuUfiUed only for much larger population sizes (indeed 
for population sizes of = 5 x 10^ we have found good agreement between the analytic 
description and the results of the simulations, over the whole range of parameters). 

We have then investigated the effects on the fluctuation power spectrum of relaxing 2 or 3 
in order to account for more realistic contact networks of the population or disease recovery 
proflles. These are examples spatial and temporal correlations that are always present in 
real interacting systems. In either case, the approximate analytic description is no longer 
valid, and we must resort to systematic simulations. However, we have found fluctuation 
power spectra that are dominantly resonant like, which suggests that the basic mechanism 
at play is still resonance with demographic stochasticity. 

Instead of assuming 2, we have modelled the mixing network of the populations as a 
dynamic small world network, and considered several different values of the small world 
parameter p that measures the degree of randomness of the contact network. We have 
found that, as p decreases, the resonant fluctuations become larger and more coherent over 
the whole range of parameters that are relevant for the description of childhood infections, 
and the dominant frequency of these fluctuations is shifted towards smaller frequencies. 

Instead of 3, we have considered the opposite (strongly correlated) limit of deterministic 
recovery at the end of the infectious period. Clearly, realistic recovery proflles will lie 
between these two extremes. We have found, again, resonant fluctuations that are much 
larger, and more coherent than when recovery occurs randomly at a constant rate, and that 
the dominant frequency is shifted to larger frequencies. 

The main conclusion is then that the importance of the phenomenon reported in [l^ to 
the description of the long term dynamics of childhood diseases is enhanced when the model 
is modifled to include more realistic assumptions (correlations), either on the populations 
contact patterns (spatial correlations) or on the disease recovery proflle (time correlations). 
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Our results apply to a large region in the epidemiological and demographic parameter space, 
and in this sense they are applicable in general to the modelling of infectious disease dynam- 
ics. They show that stochasticity may play the leading role in determining the incidence 
temporal patterns, through resonance of internal noise with the dynamics that governs the 
system in the limit of an infinite population. 

Our results also show that the analytic theory developed in provides an overall good 
description for well mixed populations and diseases with gradually decaying recovery profiles, 
and a useful guideline for the case of populations with non trivial contact networks. When 
these contact networks include a large fraction, of 50% or more, of random connections, 
an effective theory based on the analytic treatment of the randomly mixed case with a 
correction of the transmissibility to account for infective screening gives a good description 
of the incidence fluctuations spectrum. In highly correlated contact networks, however, the 
overall amplitude of the fluctuations is much larger than that predicted by this effective 
theory. 

The qualitative picture for the dependence of demographic stochasticity on the system 
parameters that emerges from our results is more subtle than what one would expect from 
linear perturbation analysis, either of the deterministic model ([1]) or even of the full stochas- 
tic description First, in the presence of correlations, the dominant frequency of these 
resonant organized fluctuations differs significantly from the natural frequency of the de- 
terministic description. The fact that this naive expectation may seem to lead to a good 
description, as argued in [g^, must be attributed to the cancelation of several missing effects. 
Second, the amplitude of these fluctuations decreases with A^, or more precisely with -\/iV, 
but finite size effects will be important for realistic population sizes, of the order of 10^, 
except for diseases with an epidemiological profile similar to that of pertussis. Third, the 
amplitude of the fluctuations increases when Rq decreases, and, for fixed Rq, it decreases as 
7 decreases. This additional dependence of the amphtude of the fluctuations on the average 
infectious period is unaccounted for by the analytic description (JSj). 

On more conceptual grounds, the finding that in finite, discrete populations internal noise 
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together with correlations produces sustained incidence oscillations of significant amplitude 
all over the parameter region that includes childhood infectious diseases is of importance for 
the long-standing controversy in epidemiology and ecology as to the driving mechanisms of 



the pervasive noisy oscillations observed in these systems [46|]. Whether these are mainly 
intrinsic or external seems to depend not only on the model's nonlinearities but also on the 
correlations between the systems's units, which most traditional approaches neglect. 

As the need for spatial models of infectious disease transmission is increasingly acknowl- 
edged, there are different modelling strategies that try to reconcile explicit spatial represen- 
tation with computational costs and the informatin available on the patterns of contact of 
the populations (see and refences therein). In the context of this discussion, it should 
be stressed that the kind of intrinsic stochastic effects highlighted in the present paper is 
specific of models that explicitely represent individuals (or small units such as households), 
and that computationally lighter, coarse-grained models will miss this phenomenology. 

V. APPENDIX 

All simulations were carried out over a population of individuals arranged on a square 
lattice with periodic boundaries where each site connects to its 12 closest neighbours. A 
random fraction p of these local links are replaced at each step with random long-distance 
links. The population size was 10^ individuals for all simulations, with the exception of the 
results of Figure 4 where a population of 50 x 10^ individuals was also considered. 

The simulations implemented the stochastic process of the SIR model, or its modification 
to include deterministic, instead of stochastic, recovery, on this network. For the case of 
stochastic recovery, we have used the efficient algorithm for stochastic processes in spatially 



structured systems described in 48|]. Local and long range infections are dealt with sepa- 
rately, and a local (resp. long range) infection event occurs with probability 1 — p (resp. p). 
For p = 1, the presence of spatial structure is irrelevant, and the algorithm reduces to the 
application of the method of Gillespie 4^ to the stochastic process ([2]). When p < 1 and 
local infections may occur, their probability depends on the number of infected neighbours 
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of each susceptible node, and so susceptible nodes with k infected neighbours, k — 0, 12 
are treated as separate classes. For the case of stochastic recovery, this algorithm is ex- 
tended with the inclusion of a linked list of infected nodes. Each entry in the list records the 
position of an infected node in the lattice, along with the timestep number (counting from 
at the start of the run) at which this particular node should recover (in case it has not 
been removed previously by a stochastic death transition). The list is ordered by recovery 
timestep number; that is, it always contains the oldest nodes at the start. In this way, 
recovery is efficiently performed by checking, at the beginning of each timestep, which nodes 
are due for recovery and replacing them in the lattice with recovered nodes 

The numerical results of Figures 3, 5 and 8 were obtained through systematic simulations 
of the SIR model for a set of 21 x 21 evenly spaced points covering the {/i/j, P/j) plane, for 
the region depicted in the figures (and keeping = 5.5 x 10"^ fixed). Four separate sets of 
simulations were carried out, one for each of the specific cases considered in the main text, 
namely, the SIR model with stochastic recovery for p = 1.0 and for p = 0.6, the SIR model 
with stochastic recovery for (/i/7,/9/7) = (0.00116,17.0) and several values of p, and the 
SIR model with deterministic recovery and p = 1.0. 

For each sample point, 500 independent runs of the model were taken, each lasting for 
50000 timesteps or days, of which only the last 2^^ are used; approximately 17000 timesteps 
are discarded at the begining, to allow each run to 'settle down' to its steady state. If less 
than 200 of these 500 runs survive (that is, finish with more than zero infectives), then 
this particular sample point is discarded (shown as grey in the plots); otherwise, the power 
spectral densities of the surviving runs are computed (using an FFT routine) and averaged 
together. The remaining power spectrum plots shown in the Figures were obtained in the 
same way, except those of Figure 4. (a), in which 2^^ timesteps long stationary timeseries 
were used, 8 times longer than those used for other plots, to provide the increased frequency 
resolution necessary to check the perfect agreement with the the analytical prediction (3) 
for 50 X 10^ individuals. 

From the final averaged spectrum, the quantities A, c and s are computed. The overall 
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amplification, which is defined ets A — ^ J^°° Pi{co) dcu, where Piiou) is the power spectrum 
density of the scaled fluctuations, is numerically approximated by a sum over the squared 
moduli of the scaled coefficients Z^/VN of the discrete Fourier transform of the timeseries. 
Since the sampling interval is 1 day, the band width of the signal is a; e [— tt, tt] , and we 
have 

1 f+'^ 1 1 -^"^ 

A^- P(uj)(hj = - P(Trk/L)dk^ -Y^^-^ 
TT Jo L Jo L N 

where k is the integer index of the Fourier coefficient and L the sample length (here L — 2^^). 
The coherence is defineded as c = Ap/A, where Ap is the integral of the power spectrum 
density of the scaled fluctuations over the frequency range that corresponds to ±10% of the 
period of the dominant peak. The numerical value of Ap is obtained in the same manner 
as A from the squared moduli of the scaled coefficients Zk/y/N, replacing the summation 
limits appropriately. 

The position of the dominant peak Upeak is found from the largest Fourier coefficient, and 
is also used to compute the shift factor s — {cOpeak ~ VD) / VD shown in Figure 7. 
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